!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!
!      This is a simple Python wrapper for the shtools library written by Mark
!      Wieczorek (http://shtools.ipgp.fr/). It uses the f2py compiler to create
!      the shtools.so shared object that can be opened from python. I tried to
!      make the wrapper as "pythonic" as possible, i.e. get the array shapes
!      etc. automatically from input arrays. This wrapper is far from complete.
!      After installation and setup of "f2py" (usually comes with
!      numpy/scipy)(http://www.scipy.org/F2py) compile with: f2py
!      -I<SHTOOLS_INCLUDE_DIR> -L<SHTOOLS_LIB_DIR> -lSHTOOLS2.8 \ -lfftw3 -lm
!      --f90flags="-m64 -fPIC" --f77flags="-m64 -fPIC" -c -m shtools wrapper.f90
!
!      in python type: >>> import shtools. Examine content with shtools.<TAB>
!      and use ipythons ? operator with shtools.<subroutine>? . You can also use
!      pythons dir(shtools) as well as print shtools.pymakegriddh.__doc__ to get
!      info on input and output variables. Also check out SHTOOLS documentary
!      for more information. We use geodesy + no Condon-Shortley Phase. To
!      change it, you have to change the subroutine below and recompile.
!
!      by Matthias Meschede (2013)
!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

subroutine pymakegriddh(coeffs,nlat,lmax_coeffs,grid)
    use shtools
    implicit none
    double precision :: interval
    integer :: n, nlat, lmax_grid, lmax_coeffs
    double precision grid(nlat,2*nlat)
    double precision coeffs(2,lmax_coeffs+1,lmax_coeffs+1)
    !f2py intent(out) grid
    !f2py intent(in) coeffs
    !f2py integer intent(in),depend(coeffs),optional :: lmax_coeffs = shape(coeffs,1)-1
    !f2py integer intent(in),depend(lmax_coeffs),optional :: nlat = 2*(lmax_coeffs +1)

    lmax_grid = nlat/2-1

    call MakeGridDH(grid,n,coeffs,lmax_grid,lmax_calc=lmax_coeffs,sampling=2)
end subroutine

subroutine pyshmultiply(shout, sh1, lmax1, sh2, lmax2)
    use shtools
    implicit none
    double precision shout(2,lmax1+lmax2+1, lmax1+lmax2+1)
    double precision sh1(2,lmax1+1,lmax1+1)
    double precision sh2(2,lmax2+1,lmax2+1)
    integer lmax1,lmax2
    !f2py intent(out) shout
    !f2py intent(in) sh1, sh2
    !f2py integer intent(in), depend(sh1), optional :: lmax1 = shape(sh1,1)-1
    !f2py integer intent(in), depend(sh2), optional :: lmax2 = shape(sh2,1)-1
    call SHMultiply(shout, sh1, lmax1, sh2, lmax2, 0)
end subroutine

subroutine pyshmultitaper_direct(mtse, mtse_dev, coeffs, lmax_coeffs, lmaxt, k, theta0, clat, clon)
    !this function performs a multitaper localized spectral analysis on a
    !cartesian box
    use shtools
    implicit none
    integer lmax_coeffs, lmaxt, k
    double precision coeffs(2,lmax_coeffs+1,lmax_coeffs+1)
    double precision mtse(lmax_coeffs - lmaxt + 1)
    double precision theta0, clat, clon
    !local variables
    double precision PI
    double precision tapers(lmaxt+1, (lmaxt+1)**2)
    double precision eigenvalues( (lmaxt+1)**2 )
    integer taper_order( (lmaxt+1)**2 )
    double precision mtse_dev(lmax_coeffs - lmaxt + 1)
    !f2py intent(out) mtse,mtse_dev
    !f2py intent(in) lmaxt,k
    !f2py intent(in) clat,clon
    !f2py integer intent(in),depend(coeffs),optional :: lmax_coeffs = shape(coeffs,1)-1
    call SHReturnTapers(theta0, lmaxt, tapers, eigenvalues, taper_order)
    call SHMultiTaperSE(mtse, mtse_dev, coeffs, lmax_coeffs, tapers, taper_order, &
                                 lmaxt, k, lat = clat, lon = clon)
end subroutine

subroutine pyshmultitaperse(mtse, mtse_dev, coeffs, lmax_coeffs, tapers, taper_order, &
                                lmaxt, ntapers, k, clat, clon)
    use shtools
    implicit none
    integer lmax_coeffs, lmaxt, k, ntapers
    double precision coeffs(2,lmax_coeffs+1,lmax_coeffs+1)
    double precision mtse(lmax_coeffs - lmaxt + 1)
    double precision clat, clon
    double precision tapers(lmaxt+1, ntapers)
    integer taper_order( ntapers )
    double precision mtse_dev(lmax_coeffs - lmaxt + 1)
    !f2py intent(out) mtse,mtse_dev
    !f2py integer intent(in),depend(tapers),optional :: lmaxt = shape(tapers,0)-1
    !f2py integer intent(in),depend(tapers),optional :: ntapers = shape(tapers,1)
    !f2py intent(in) k
    !f2py integer intent(in),depend(coeffs),optional :: lmax_coeffs = shape(coeffs,1)-1
    !f2py intent(in) tapers
    !f2py intent(in) taper_order
    !f2py intent(in) clat,clon
    call SHMultiTaperSE(mtse, mtse_dev, coeffs, lmax_coeffs, tapers, taper_order, lmaxt, k, &
                                                                       lat = clat, lon = clon)
end subroutine

subroutine pyshmultitapercse(mtse, mtse_dev, coeffs1, coeffs2, lmax_coeffs, &
                             tapers, taper_order, lmaxt, ntapers, k, theta0, clat, clon)
    use shtools
    implicit none
    integer lmax_coeffs, lmaxt, k, ntapers
    double precision coeffs1(2,lmax_coeffs+1,lmax_coeffs+1)
    double precision coeffs2(2,lmax_coeffs+1,lmax_coeffs+1)
    double precision mtse(lmax_coeffs - lmaxt + 1)
    double precision theta0, clat, clon
    double precision tapers(lmaxt+1, ntapers)
    integer taper_order(ntapers)
    double precision mtse_dev(lmax_coeffs - lmaxt + 1)
    !f2py intent(out) mtse,mtse_dev
    !f2py integer intent(in),depend(tapers),optional :: lmaxt = shape(tapers,0)-1
    !f2py integer intent(in),depend(tapers),optional :: ntapers = shape(tapers,1)
    !f2py intent(in) k
    !f2py intent(in) coeffs1, coeffs2
    !f2py integer intent(in),depend(coeffs1),optional :: lmax_coeffs = shape(coeffs1,1)-1
    !f2py intent(in) tapers
    !f2py intent(in) taper_order
    !f2py intent(in) theta0,clat,clon
    call SHMultiTaperCSE(mtse, mtse_dev, coeffs1, lmax_coeffs, coeffs2, lmax_coeffs, &
                                   tapers, taper_order, lmaxt, k, lat = clat, lon = clon)
end subroutine

subroutine pyshreturntapersmap(tapers, eigenvalues, dh_mask, n, lmaxt, ntapers)
    use shtools
    implicit none
    double precision tapers((lmaxt+1)*(lmaxt+1), ntapers)
    double precision eigenvalues(ntapers)
    integer lmaxt,n,ntapers
    integer dh_mask(N,2*N)
    !f2py intent(out) tapers, eigenvalues
    !f2py intent(in) dh_mask, n, lmaxt, ntapers
    !f2py integer intent(in), depend(dh_mask), optional :: n = shape(dh_mask,0)
    call SHReturnTapersMap(tapers, eigenvalues, dh_mask, n, 2, lmaxt, ntapers)
end subroutine

subroutine pyshlocalizedadmitcorr(tapers, taper_order, lwin, clat, clon,&
    coeffs1, coeffs2, lmax,admit,corr, k, admit_error, corr_error)
    use shtools
    implicit none
    integer lwin, lmax, k, i, j 
    integer taper_order(k)
    double precision tapers(lwin+1, k)
    double precision coeffs1(2,lmax+1,lmax+1)
    double precision coeffs2(2,lmax+1,lmax+1)
    double precision clat,clon
    double precision, dimension(lmax-lwin+1) :: admit,corr,admit_error,corr_error
    !f2py intent(in) tapers, taper_order
    !f2py intent(in) coeffs1, coeffs2
    !f2py intent(out) admit, corr, admit_error, corr_error
    !f2py integer intent(in),depend(tapers),optional :: k = shape(tapers,1)
    !f2py integer intent(in),depend(tapers),optional :: lwin = shape(tapers,0)-1
    !f2py integer intent(in),depend(coeffs1),optional :: lmax = shape(coeffs1,1)-1
    print *, 'localized admittence and correlation function got shapes:'
    print *, 'coeffs1:', shape(coeffs1)
    print *, 'coeffs2:', shape(coeffs2)
    call SHLocalizedAdmitCorr(tapers, taper_order, lwin, clat, clon, coeffs1,&
                             coeffs2, lmax, admit, corr, k, admit_error, corr_error)
end subroutine

subroutine pyshmtvaropt(l, tapers, taper_order, lmaxt, k, sff, var_opt, var_unit, weight_opt, &
                                                        unweighted_covar, nocross)
    use shtools
    implicit none
    integer l,k
    double precision tapers(lmaxt+1,k)
    double precision sff()
end subroutine

subroutine pyshvectortocilm(vector, cilm, lmax)
    use shtools
    implicit none
    double precision vector((lmax+1)*(lmax+1))
    double precision cilm(2,lmax+1,lmax+1)
    integer lmax
    !f2py intent(in)  lmax,vector
    !f2py intent(out) cilm
    call SHVectorToCilm(vector, cilm, lmax)
end subroutine

subroutine pyshreturntapers(theta0, lmax, tapers, eigenvalues, taper_order)
    use shtools
    implicit none
    integer lmax
    double precision theta0
    double precision tapers(lmax+1, (lmax+1)*(lmax+1))
    double precision eigenvalues((lmax+1)*(lmax+1))
    integer taper_order( (lmax+1)*(lmax+1) )
    !f2py intent(in) theta0, lmax
    !f2py intent(out) tapers, eigenvalues, taper_order
    call SHReturnTapers(theta0, lmax, tapers, eigenvalues, taper_order)
end subroutine

subroutine pyshbiask(tapers, lmaxt, k, k_tap, inspec, lmax_in, outspec)
    use shtools
    implicit none
    double precision tapers(lmaxt+1, k_tap)
    double precision inspec(lmax_in+1)
    double precision outspec(lmax_in+1)
    integer lmaxt, k_tap, k, lmax_in
    !f2py intent(in) tapers, inspec, k
    !f2py integer intent(in), depend(tapers), optional :: lmaxt   = shape(tapers,0)-1
    !f2py integer intent(in), depend(tapers), optional :: k_tap   = shape(tapers,1)
    !f2py integer intent(in), depend(inspec), optional :: lmax_in = shape(inspec,0)-1
    !f2py intent(out) outspec
    call SHBiasK(tapers, lmaxt, k, inspec, lmax_in, outspec)
end subroutine

subroutine pyshmtdebias(mtdebias, mtspectra, lmax, nl, tapers, lwin, k, lmid,n )
    use shtools
    implicit none
    double precision mtdebias(2, (lmax+2)/nl )
    double precision mtspectra(2, lmax+1)
    double precision tapers(lwin+1, k)
    double precision lmid( (lmax+2)/nl )
    integer lmax, nl, lwin, k, n
    !f2py intent(in) mtspectra
    !f2py intent(in) nl
    !f2py intent(in) tapers
    !f2py integer intent(in), depend(mtspectra), optional :: lmax = shape(mtspectra,1)-1
    !f2py integer intent(in), depend(tapers), optional :: lwin = shape(tapers,0)-1
    !f2py integer intent(in), depend(tapers), optional :: k    = shape(tapers,1)
    !f2py intent(out) mtdebias
    !f2py intent(out) lmid
    !f2py intent(out) n
    call SHMTDebias(mtdebias, mtspectra, lmax, tapers, lwin, k, nl, lmid, n)
end subroutine

subroutine pyshexpanddh(grid,nlat,coeffs,lmax_coeffs)
    use shtools
    implicit none
    double precision grid(nlat,2*nlat)
    double precision coeffs(2,lmax_coeffs+1,lmax_coeffs+1)
    integer nlat,lmax_coeffs,lmax_grid
    !f2py intent(out) coeffs
    !f2py intent(in) grid
    !f2py integer intent(hide), optional, depend(grid) :: nlat = shape(grid,0)
    !f2py integer intent(in), optional, depend(nlat):: lmax_coeffs = nlat/2-1

    call SHExpandDH(grid,nlat,coeffs,lmax_grid,lmax_calc=lmax_coeffs,sampling=2)
end subroutine

subroutine pymakegridpoint(coeffs,lmax_calc,lmax_coeffs,lat,lon,return_value)
    use shtools
    implicit none
    integer lmax_calc,lmax_coeffs
    double precision coeffs(2,lmax_coeffs+1,lmax_coeffs+1)
    double precision lat, lon
    double precision return_value
    !f2py intent(in) coeffs
    !f2py intent(in) lat
    !f2py intent(in) lon
    !f2py intent(out) return_value
    !f2py integer intent(in),depend(coeffs),optional :: lmax_coeffs = shape(coeffs,1)-1
    !f2py integer intent(in),depend(lmax_coeffs),optional :: lmax_calc = lmax_coeffs
    return_value = MakeGridPoint(coeffs, lmax_calc, lat, lon)
end subroutine

subroutine pyshpowerspectrum(coeffs,lmax,spectrum)
    use shtools
    implicit none
    double precision coeffs(2,lmax+1,lmax+1)
    double precision spectrum(lmax+1)
    integer lmax
    !f2py intent(out) spectrum
    !f2py intent(in)  coeffs
    !f2py integer intent(hide), depend(coeffs) :: lmax = shape(coeffs,1) - 1

    call SHPowerSpectrum(coeffs, lmax, spectrum)
end subroutine

subroutine pyshcrosspowerspectrum(coeffs1, coeffs2, lmax, cspectrum)
    use shtools
    implicit none
    double precision coeffs1(2,lmax+1,lmax+1)
    double precision coeffs2(2,lmax+1,lmax+1)
    double precision cspectrum(lmax+1)
    integer lmax
    !f2py intent(in) coeffs1
    !f2py intent(in) coeffs2
    !f2py integer intent(hide),depend(coeffs1) :: lmax = shape(coeffs1,1) - 1
    !f2py intent(out) cspectrum

    call SHCrossPowerSpectrum(coeffs1,coeffs2,lmax,cspectrum)
end subroutine

subroutine pyshcrosspowerspectrumdensity(coeffs1, coeffs2, lmax, cspectrum)
    use shtools
    implicit none
    double precision coeffs1(2,lmax+1,lmax+1)
    double precision coeffs2(2,lmax+1,lmax+1)
    double precision cspectrum(lmax+1)
    integer lmax
    !f2py intent(in) coeffs1
    !f2py intent(in) coeffs2
    !f2py integer intent(hide),depend(coeffs1) :: lmax = shape(coeffs1,1) - 1
    !f2py intent(out) cspectrum

    call SHCrossPowerSpectrumDensity(coeffs1,coeffs2,lmax,cspectrum)
end subroutine

subroutine pyshadmitcorr(coeffs1,coeffs2,lmax,admit,corr, admit_error)
    use shtools
    implicit none
    double precision coeffs1(2,lmax+1,lmax+1)
    double precision coeffs2(2,lmax+1,lmax+1)
    double precision admit(lmax+1)
    double precision corr(lmax+1)
    double precision admit_error(lmax+1)
    integer lmax
    !f2py intent(in) coeffs1
    !f2py intent(in) coeffs2
    !f2py integer intent(hide),depend(coeffs1) :: lmax = shape(coeffs1,1) - 1
    !f2py intent(out) admit
    !f2py intent(out) corr
    !f2py intent(out) admit_error

    call SHAdmitCorr(coeffs1, coeffs2, lmax, admit, corr, admit_error=admit_error)
end subroutine

subroutine pyshconfidence(l, corr)
    implicit none
    integer l
    double precision corr
    !f2py intent(in) l
    !f2py intent(in) corr

    call SHConfidence(l,corr)
end subroutine

subroutine pyshrotaterealcoef(coeffs_rot, coeffs, lmax, angles)
    use shtools
    implicit none
    !---- external ----
    double precision coeffs_rot(2,lmax+1,lmax+1)
    double precision coeffs(2,lmax+1,lmax+1)
    double precision angles(3)
    double precision DJ(lmax+1,lmax+1,lmax+1)
    integer lmax
    !f2py intent(in) coeffs
    !f2py intent(in) angles
    !f2py integer intent(hide),depend(coeffs) :: lmax = shape(coeffs,1) - 1
    !f2py intent(out) coeffs_rot

    call djpi2(DJ,lmax)
    call SHRotateRealCoef(coeffs_rot, coeffs, lmax, angles, DJ)
end subroutine

subroutine pyplmbar(p, lmax, cost)
    use shtools
    implicit none
    integer lmax
    double precision :: cost
    double precision :: p( (lmax+1)*(lmax+2)/2 ) !SHtools array (contains all ls&ms)
    double precision :: value
    !f2py intent(in) lmax,cost
    !f2py intent(out) p

    call PlmBar ( p, lmax, cost )
end subroutine

subroutine pyplmindex(l,m,ind)
    use shtools
    implicit none
    integer l,m,ind
    !f2py intent(in) l,m
    !f2py intent(out) ind
    ind = PlmIndex(l,m)
end subroutine

